#'
#'  Project title:     'Sovereign Risk and Government Change: Elections, Ideology and Experience'
#'  Authors:           Sarah M. Brooks; Raphael Cunha; Layna Mosley
#'  File description:  Heteroskedastic regression analysis of monthly EMBI and CDS spreads
#'  Output:            Tables 4-5; Figures 1-2; + robustness checks
#'

# my_packages <- c("tidyverse", "stargazer",
#                  "gtable", "gridExtra", "rio",
#                  "zoo", "reshape2", "grid",
#                  "gridExtra", "countrycode")
# 
# install.packages(my_packages)

library(tidyverse)
library(utils)
library(rio)
library(reshape2)
library(gtable)
library(stargazer)
library(gridExtra)
library(zoo)
library(countrycode)

# Set directory

SAVEDIR <- "~/.../"

# Load EMBIG and CDS Spreads Dataset

load(file.path(SAVEDIR, "Data/Monthly", "EMBIG and CDS Dataset.RData"))

# Summary statistics

sum_stat <- spread.dat %>%
  ungroup() %>%
  dplyr::select(embispread, d_embispread, cdsspread, d_cdsspread,
                monthsoffc, execleft, election_window,
                cab_annual_gdp, extdebt_gni_annual, growth_annual, shortdebt_reserves_annual, inflation, d_inflation,
                growth_qtr, cab_qtr_gdp, extdebt_gdp_qtr,
                treasury10y, d_treasury10y, vix, d_vix, energy_index, d_energy_index, eqprem, d_eqprem, globaldefault,
                embispreadREGIONdiff, d_embispreadREGIONdiff, cdsspreadREGIONdiff, d_cdsspreadREGIONdiff,
                kaopen,
                fxrate, d_fxrate,
                left_region, baker_ideol,
                SPratingnumeric, d_SPratingnumeric, moodysratingnumeric, d_moodysratingnumeric, fitchratingnumeric, d_fitchratingnumeric,
                v2x_libdem, polconv) %>%
  as.data.frame() %>%
  stargazer(digits = 2)


# Estimate heteroskedastic regression models ------------------------------

# Load 'hetreg' and 'woolARtest' functions

source(file.path(SAVEDIR, "2b_functions.R"))

# Add interactions for right governments

spread.dat$execrightXmonthsoffc <- spread.dat$execright * spread.dat$monthsoffc
spread.dat$execrightXloss <- spread.dat$execright * spread.dat$loss

# Lag DVs

spread.dat <- spread.dat %>%
  group_by(iso3c) %>%
  mutate(l_d_embispread = dplyr::lag(d_embispread),
         l_d_cdsspread = dplyr::lag(d_cdsspread))


# Main results (Table 4) --------------------------------------------------

fe.var <- "country" # fixed-effect unit
y.var <- "d_embispread" # outcome measure

x.var <- c("l_d_embispread",
           "monthsoffc", "execleft")
x.het <- c("monthsoffc", "execleft")
mod.1 <- hetreg(y.var, x.var, x.int = NULL, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("l_d_embispread",
           "monthsoffc")
x.int <- c("monthsoffc", "execleft")
x.het <- c("monthsoffc")
mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("l_d_embispread",
           "election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("election_window", "d_vix", "kaopen")
mod.3 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("l_d_embispread",
           "election_window", "execright", "execrightXmonthsoffc",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("election_window", "d_vix", "kaopen", "execright", "execrightXmonthsoffc")
mod.4 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

y.var <- "d_cdsspread" # Outcome measure

x.var <- c("monthsoffc", "execleft")
x.het <- c("monthsoffc", "execleft")
mod.5 <- hetreg(y.var, x.var, x.int = NULL, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("monthsoffc")
x.int <- c("monthsoffc", "execleft")
x.het <- c("monthsoffc")
mod.6 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("election_window", "d_vix", "kaopen")
mod.7 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window", "execright", "execrightXmonthsoffc",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("election_window", "d_vix", "kaopen", "execright", "execrightXmonthsoffc")
mod.8 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

models <- list(mod.1, mod.2, mod.3, mod.4,
               mod.5, mod.6, mod.7, mod.8)

models_table <- map(models, function (x) {x[["cl.SE.table"]]})
n_obs <- map(models, function (x) {x[["nobs"]]})
n_obs <- c("Observations", do.call("cbind", n_obs))
n_countries <- map(models, function (x) {x[["countries"]]})
n_countries <- c("Countries", do.call("cbind", n_countries))
aic <- map(models, function (x) {x[["aic"]]})
aic <- c("AIC", do.call("cbind", aic))
country_fe <- c("Country Fixed Effects", rep("\\checkmark", length(models)))

wool_test <- map(models, function (x) {
  
  woolARtest(model_residuals = residuals(x[["mod"]]),
             model_data = x[["data"]],
             id_var = "country",
             time_var = "month_year")
  })

AR_test1 <- map(wool_test, function (x) {format(x$statistic, digits = 1, nsmall = 2, scientific = FALSE)})
AR_test1 <- c("Wooldridge AR(1) Test $F$-Stat.", do.call("cbind", AR_test1))
AR_test2 <- map(wool_test, function (x) {paste("(\\textit{p} = ", format(x$p.value, digits = 1, nsmall = 2, scientific = FALSE), ")", sep = "")})
AR_test2 <- c("", do.call("cbind", AR_test2))

stargazer(models_table,
          add.lines = list(n_obs, n_countries, country_fe, aic,
                           AR_test1, AR_test2),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)

# Plot predicted volatility under left governments by months in office

# Function for plotting interaction in volatility
# equation using simulated Clarify CIs

clariplot <- function(mod){
  
  # Draw parameters from multivariate distribution
  S <- 10000
  set.seed(98765)
  
  sim.coef <- mvrnorm(n = S, mu = mod[["betas"]][,"coefficients"], Sigma = mod[["sandwich"]])
  sim.sigma.coef <- as.matrix(dplyr::select(as.data.frame(sim.coef), starts_with("sigma_")))
  
  # Set covariates at desired values
  monthsoffc <- seq(1, 96, by = 1)
  execleft.1 <- 1
  execleft.0 <- 0
  vix <- median(unique(spread.dat[, c("month_year", "d_vix")])$d_vix)
  kaopen <- median(spread.dat$kaopen, na.rm = TRUE)
  election_window <- 0

  out.1 <- matrix(NA, nrow = length(monthsoffc), ncol = 5)
  out.0 <- matrix(NA, nrow = length(monthsoffc), ncol = 5)

  for(i in 1:length(monthsoffc)){
    
    x.1 <- cbind(1, monthsoffc[i], execleft.1, vix, kaopen, election_window,
                 monthsoffc[i] * execleft.1)
      
    pred.sim.1 <- exp(sim.sigma.coef %*% t(x.1))
      
    x.0 <- cbind(1, monthsoffc[i], execleft.0, vix, kaopen, election_window,
                 monthsoffc[i] * execleft.0)
      
    pred.sim.0 <- exp(sim.sigma.coef %*% t(x.0))
      
    pred.sim <- pred.sim.1 - pred.sim.0
      
    out.1[i,] <- c(1, i, mean(pred.sim.1, na.rm = TRUE), quantile(pred.sim.1, c(0.025, 0.975), na.rm = TRUE))
    out.0[i,] <- c(0, i, mean(pred.sim.0, na.rm = TRUE), quantile(pred.sim.0, c(0.025, 0.975), na.rm = TRUE))
    
    rm(x.1, x.0, pred.sim.1, pred.sim.0, pred.sim)
    }

  out <- as.data.frame(rbind(out.1, out.0))
  colnames(out) <- c("left", "monthsoffc", "mean", "lb", "ub")
  
  # Get coefficient and t-stat on interaction term for annotation
  b_int <- signif(mod[["cl.SE.table"]]["sigma_monthsoffc:execleft", 1], digits = 2)
  t_int <- format(round(mod[["cl.SE.table"]]["sigma_monthsoffc:execleft", 3], digits = 2), nsmall = 2)

  return(list(out = out, b_int = b_int, t_int = t_int))
}

# Common plot theme

p_theme <- theme_bw(base_family = "Myriad Pro") +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        plot.title =  element_text(hjust = 0.5),
        plot.subtitle =  element_text(hjust = 0.5),
        #panel.grid.minor = element_blank(),
        panel.grid = element_line(linetype = 3),
        strip.background = element_blank(),
        strip.text = element_text(face = "bold", size = 10))

plot_dat <- clariplot(mod.3)

up_limit <- max(plot_dat[["out"]]$ub)

p1 <- ggplot(plot_dat[["out"]], aes(x = monthsoffc, y = mean, group = factor(left), linetype = factor(left))) +
  geom_ribbon(aes(x = monthsoffc, ymin = lb, ymax = ub), fill = "grey80", alpha = 0.3) +
  geom_line() +
  labs(x = NULL, y = NULL) +
  scale_linetype_manual(values = c(1, 2), labels = c("Right/Center", "Left")) +
  ylim(25, up_limit) +
  ggtitle(label = "EMBI Spread",
          #subtitle = paste("Coef. on 'Left Gov. x Months in Office' is ", plot_dat[["b_int"]], "\n(t-statistic = ", plot_dat[["t_int"]], ")", sep = "")
          subtitle = bquote(atop("Coef. on" ~ italic("Left Gov.") %*% italic("Months in Office") ~ "is" ~ .(plot_dat[["b_int"]]), "t-statistic = " ~ .(plot_dat[["t_int"]])))) +
  p_theme+
  NULL

plot_dat <- clariplot(mod.7)

p2 <- ggplot(plot_dat[["out"]], aes(x = monthsoffc, y = mean, group = factor(left), linetype = factor(left))) + 
  geom_ribbon(aes(x = monthsoffc, ymin = lb, ymax = ub), fill = "grey80", alpha = 0.3) +
  geom_line() +
  labs(x = NULL, y = NULL) +
  scale_linetype_manual(values = c(1, 2), labels = c("Right/Center", "Left")) +
  ylim(25, up_limit) +
  ggtitle(label = "CDS Spread",
          #subtitle = paste("Coef. on 'Left Gov. x Months in Office' is ", plot_dat[["b_int"]], "\n(t-statistic = ", plot_dat[["t_int"]], ")", sep = ""))
          subtitle = bquote(atop("Coef. on" ~ italic("Left Gov.") %*% italic("Months in Office") ~ "is" ~ .(plot_dat[["b_int"]]), "t-statistic = " ~ .(plot_dat[["t_int"]])))) +
  p_theme +
  NULL

# Combine plots

legend <- gtable_filter(ggplot_gtable(ggplot_build(p1 + theme(legend.position = "bottom"))), "guide-box")
lheight <- sum(legend$height)
p <- arrangeGrob(p1 + theme(legend.position = "none"), p2 + theme(legend.position = "none"),
                 ncol = 2, left = "Spread volatility", bottom = textGrob("Months in office", gp = gpar(fontfamily = "Myriad Pro")))
theight <- unit(12, "points")
p <- arrangeGrob(p, legend, heights=unit.c(unit(1, "npc") - lheight*2, lheight*1.5), ncol = 1)
grid.arrange(p)


# Country and time coverage of analysis

cover_embi <- mod.1[["data"]] %>%
  mutate(date = as.yearmon(month_year, "%m-%Y")) %>%
  group_by(country) %>%
  dplyr::summarise(min_embi = as.character(min(date)),
                   max_embi = as.character(max(date))) %>%
  mutate(embi = paste(min_embi, " - ", max_embi)) %>%
  dplyr::select(country, embi) %>%
  as.data.frame()

cover_cds <- mod.4[["data"]] %>%
  mutate(date = as.yearmon(month_year, "%m-%Y")) %>%
  group_by(country) %>%
  dplyr::summarise(min_cds = as.character(min(date)),
                   max_cds = as.character(max(date))) %>%
  mutate(cds = paste(min_cds, " - ", max_cds)) %>%
  dplyr::select(country, cds) %>%
  as.data.frame()

cover <- full_join(cover_embi, cover_cds, by = "country")
cover <- cover[order(cover$country), ]

print(xtable(cover), include.rownames = FALSE) 


# Common support for 'partisanship x time in office' interaction

p <- ggplot(spread.dat[spread.dat$monthsoffc <= 120, ],
            aes(x = monthsoffc, fill = factor(execleft))) +
  geom_histogram(colour = "grey60", size = .3, alpha = 0.5, binwidth = 1) +
  labs(x = "Months in office", y = "Count") +
  theme_bw(base_family = "Myriad Pro") +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        plot.caption =  element_text(hjust = 0.5),
        panel.grid = element_blank()) +
  scale_fill_grey(labels = c("Right/Center", "Left", "Unclassified"),
                  start = 1, end = 0.4) +
  NULL
p


# Conditional on global interest rates (Table 5) --------------------------

fe.var <- "country" # fixed-effect unit

y.var <- "d_embispread" # outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft", "treasury10y")
x.het <- c("d_vix", "kaopen", "election_window")
mod.1 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)


y.var <- "d_cdsspread" # outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft", "treasury10y")
x.het <- c("d_vix", "kaopen", "election_window")
mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

models <- list(mod.1, mod.2)

models_table <- map(models, function (x) {x[["cl.SE.table"]]})
n_obs <- map(models, function (x) {x[["nobs"]]})
n_obs <- c("Observations", do.call("cbind", n_obs))
n_countries <- map(models, function (x) {x[["countries"]]})
n_countries <- c("Countries", do.call("cbind", n_countries))
aic <- map(models, function (x) {x[["aic"]]})
aic <- c("AIC", do.call("cbind", aic))
country_fe <- c("Country Fixed Effects", rep("\\checkmark", length(models)))

wool_test <- map(models, function (x) {
  
  woolARtest(model_residuals = residuals(x[["mod"]]),
             model_data = x[["data"]],
             id_var = "country",
             time_var = "month_year")
})

AR_test1 <- map(wool_test, function (x) {format(x$statistic, digits = 1, nsmall = 2, scientific = FALSE)})
AR_test1 <- c("Wooldridge AR(1) Test $F$-Stat.", do.call("cbind", AR_test1))
AR_test2 <- map(wool_test, function (x) {paste("(\\textit{p} = ", format(x$p.value, digits = 1, nsmall = 2, scientific = FALSE), ")", sep = "")})
AR_test2 <- c("", do.call("cbind", AR_test2))

stargazer(models_table,
          add.lines = list(n_obs, n_countries, country_fe, aic,
                           AR_test1, AR_test2),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)

# Plot triple interaction with interest rates

clariplot2 <- function(mod){

  S <- 10000
  set.seed(98765)
  sim.coef <- mvrnorm(n = S, mu = mod[["betas"]][,"coefficients"], Sigma = mod[["sandwich"]])
  sim.sigma.coef <- as.matrix(dplyr::select(as.data.frame(sim.coef), starts_with("sigma_")))
  
  # Set covariates at desired values
  monthsoffc <- seq(1, 96, by = 1)
  execleft.1 <- 1
  execleft.0 <- 0
  window <- 0
  vix <- median(unique(spread.dat[, c("month_year", "d_vix")])$d_vix)
  kaopen <- median(spread.dat$kaopen, na.rm = TRUE)
  election_window <- 0
  
  treasury10y <- quantile(unique(spread.dat[, c("month_year", "treasury10y")])$treasury10y, probs = c(0.25, 0.75), na.rm = TRUE)
  
  mfx.df <- list()
  
  for (t in 1:length(treasury10y)){
    
    out.1 <- matrix(NA, nrow = length(monthsoffc), ncol = 6)
    out.0 <- matrix(NA, nrow = length(monthsoffc), ncol = 6)
    
    for(i in 1:length(monthsoffc)){
      
      x.1 <- cbind(1, monthsoffc[i], execleft.1, treasury10y[t], vix, kaopen, election_window,
                   monthsoffc[i] * execleft.1, monthsoffc[i] * treasury10y[t], execleft.1 * treasury10y[t], monthsoffc[i] * execleft.1 * treasury10y[t])
      
      pred.sim.1 <- exp(sim.sigma.coef %*% t(x.1))
      
      x.0 <- cbind(1, monthsoffc[i], execleft.0, treasury10y[t], vix, kaopen, election_window,
                   monthsoffc[i] * execleft.0, monthsoffc[i] * treasury10y[t], execleft.0 * treasury10y[t], monthsoffc[i] * execleft.0 * treasury10y[t])
      
      pred.sim.0 <- exp(sim.sigma.coef %*% t(x.0))
      
      out.1[i,] <- c(treasury10y[t], 1, i, mean(pred.sim.1, na.rm = TRUE), quantile(pred.sim.1, c(0.025, 0.975), na.rm = TRUE))
      out.0[i,] <- c(treasury10y[t], 0, i, mean(pred.sim.0, na.rm = TRUE), quantile(pred.sim.0, c(0.025, 0.975), na.rm = TRUE))
      
      rm(x.1, x.0, pred.sim.1, pred.sim.0, pred.sim)
    }
    
    out <- rbind(out.1, out.0)
    mfx.df[[length(mfx.df) + 1]] <- as.data.frame(out)
    rm(out)
  }
  
  mfx.df <- do.call("rbind", mfx.df)
  colnames(mfx.df) <- c("trsy", "left", "monthsoffc", "mean", "lb", "ub")
  mfx.df$trsy <- factor(mfx.df$trsy)
  levels(mfx.df$trsy) <- c("Low U.S. interest rates", "High U.S. interest rates")
  
  # Get coefficient and t-stat on interaction term for annotation
  b_int <- signif(mod[["cl.SE.table"]]["sigma_monthsoffc:execleft:treasury10y", 1], digits = 2)                     # coef
  t_int <- format(round(mod[["cl.SE.table"]]["sigma_monthsoffc:execleft:treasury10y", 3], digits = 2), nsmall = 2)  # t-stat
  p_int <- format(round(mod[["cl.SE.table"]]["sigma_monthsoffc:execleft:treasury10y", 4], digits = 2), nsmall = 3)  # p-value

  return(list(mfx.df = mfx.df, b_int = b_int, t_int = t_int, p_int = p_int))
}

# Common plot theme

p_theme <- theme_bw(base_family = "Myriad Pro") +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(face = "bold", size = 10),
        panel.grid = element_line(linetype = 3))

plot_dat <- clariplot2(mod.1)

p1 <- ggplot(plot_dat[["mfx.df"]], aes(x = monthsoffc, y = mean, group = factor(left), linetype = factor(left))) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = "grey70", alpha = 0.3) +
  geom_line() +
  facet_wrap(~ factor(trsy)) +
  scale_linetype_manual(values = c(1, 2), labels = c("Right/Center", "Left")) +
  labs(x = NULL,
       y = "Spread volatility") +
  ggtitle(label = "EMBI Spread",
          subtitle = bquote("t-statistic on" ~ italic("Left Gov.") %*% italic("Months in Office") %*% italic("Treasury Rate") ~ "is" ~ .(plot_dat[["t_int"]]))) +
  ylim(0, NA) +
  p_theme +
  NULL

plot_dat <- clariplot2(mod.2)

p2 <- ggplot(plot_dat[["mfx.df"]], aes(x = monthsoffc, y = mean, group = factor(left), linetype = factor(left))) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = "grey70", alpha = 0.3) +
  geom_line() +
  facet_wrap(~ factor(trsy)) +
  scale_linetype_manual(values = c(1, 2), labels = c("Right/Center", "Left")) +
  labs(x = NULL,
       y = "Spread volatility") +
  ggtitle(label = "CDS Spread",
          subtitle = bquote("t-statistic on" ~ italic("Left Gov.") %*% italic("Months in Office") %*% italic("Treasury Rate") ~ "is" ~ .(plot_dat[["t_int"]]))) +
  ylim(0, NA) +
  p_theme +
  NULL

# Combine plots

legend <- gtable_filter(ggplot_gtable(ggplot_build(p1 + theme(legend.position = "bottom"))), "guide-box")
lheight <- sum(legend$height)
p <- arrangeGrob(p1 + theme(legend.position = "none"), p2 + theme(legend.position = "none"),
                 ncol = 1, bottom = textGrob("Months in office", gp = gpar(fontfamily = "Myriad Pro")))
theight <- unit(12, "points")
p <- arrangeGrob(p, legend, heights=unit.c(unit(1, "npc") - lheight*2, lheight*1.5), ncol = 1)
grid.arrange(p)


# Additional analyses and robustness checks -------------------------------

# Conditional on populism -------------------------------------------------

fe.var <- c("country") # Fixed-effect unit

y.var <- "d_embispread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "populist")
x.het <- c("d_vix", "kaopen", "election_window")

mod.1 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft", "populist")
x.het <- c("d_vix", "kaopen", "election_window")

mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)


y.var <- "d_cdsspread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "populist")
x.het <- c("d_vix", "kaopen", "election_window")

mod.3 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft", "populist")
x.het <- c("d_vix", "kaopen", "election_window")

mod.4 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

stargazer(mod.1[["cl.SE.table"]], mod.2[["cl.SE.table"]], mod.3[["cl.SE.table"]], mod.4[["cl.SE.table"]],
          add.lines = list(c("Observations", mod.1[["nobs"]], mod.2[["nobs"]], mod.3[["nobs"]], mod.4[["nobs"]]),
                           c("Countries", mod.1[["countries"]], mod.2[["countries"]], mod.3[["countries"]], mod.4[["countries"]]),
                           c("Country Fixed Effects", "Yes", "Yes",  "Yes", "Yes"),
                           c("AIC", mod.1[["aic"]], mod.2[["aic"]], mod.3[["aic"]], mod.4[["aic"]])),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)

# Plot triple interaction of partisanship, populism, and time in office

clariplot3 <- function(mod){
  
  S <- 10000
  set.seed(98765)
  sim.coef <- mvrnorm(n = S, mu = mod[["betas"]][,"coefficients"], Sigma = mod[["sandwich"]])
  sim.sigma.coef <- as.matrix(dplyr::select(as.data.frame(sim.coef), starts_with("sigma_")))
  
  # Set covariates at desired values
  monthsoffc <- seq(1, 96, by = 1)
  execleft.1 <- 1
  execleft.0 <- 0
  populist <- c(0, 1)
  window <- 0
  vix <- median(unique(spread.dat[, c("month_year", "d_vix")])$d_vix)
  kaopen <- median(spread.dat$kaopen, na.rm = TRUE)
  election_window <- 0
  
  mfx.df <- list()
  
  for (p in 1:length(populist)){
    
    out.1 <- matrix(NA, nrow = length(monthsoffc), ncol = 6)
    out.0 <- matrix(NA, nrow = length(monthsoffc), ncol = 6)
    
    for(i in 1:length(monthsoffc)){
      
      x.1 <- cbind(1, monthsoffc[i], execleft.1, populist[p], vix, kaopen, election_window,
                   monthsoffc[i] * execleft.1, monthsoffc[i] * populist[p], execleft.1 * populist[p], monthsoffc[i] * execleft.1 * populist[p])
      
      pred.sim.1 <- exp(sim.sigma.coef %*% t(x.1))
      
      x.0 <- cbind(1, monthsoffc[i], execleft.0, populist[p], vix, kaopen, election_window,
                   monthsoffc[i] * execleft.0, monthsoffc[i] * populist[p], execleft.0 * populist[p], monthsoffc[i] * execleft.0 * populist[p])
      
      pred.sim.0 <- exp(sim.sigma.coef %*% t(x.0))
      
      out.1[i,] <- c(populist[p], 1, i, mean(pred.sim.1, na.rm = TRUE), quantile(pred.sim.1, c(0.025, 0.975), na.rm = TRUE))
      out.0[i,] <- c(populist[p], 0, i, mean(pred.sim.0, na.rm = TRUE), quantile(pred.sim.0, c(0.025, 0.975), na.rm = TRUE))
      
      rm(x.1, x.0, pred.sim.1, pred.sim.0, pred.sim)
    }
    
    out <- rbind(out.1, out.0)
    mfx.df[[length(mfx.df) + 1]] <- as.data.frame(out)
    rm(out)
  }
  
  mfx.df <- do.call("rbind", mfx.df)
  colnames(mfx.df) <- c("populist", "left", "monthsoffc", "mean", "lb", "ub")
  mfx.df$populist <- factor(mfx.df$populist)
  levels(mfx.df$populist) <- c("Non-populist", "Populist")
  
  # Get coefficient and t-stat on interaction term for annotation
  b_int <- signif(mod[["cl.SE.table"]]["sigma_monthsoffc:execleft:populistPopulist", 1], digits = 2)                     # coef
  t_int <- format(round(mod[["cl.SE.table"]]["sigma_monthsoffc:execleft:populistPopulist", 3], digits = 2), nsmall = 2)  # t-stat
  p_int <- format(round(mod[["cl.SE.table"]]["sigma_monthsoffc:execleft:populistPopulist", 4], digits = 2), nsmall = 3)  # p-value
  
  return(list(mfx.df = mfx.df, b_int = b_int, t_int = t_int, p_int = p_int))
}

# Common theme for plots

p_theme <- theme_bw(base_family = "Myriad Pro") +
  theme(legend.position = "bottom",
        plot.caption =  element_text(hjust = 0.5),
        legend.title = element_blank(),
        strip.background = element_blank(),
        strip.text = element_text(face = "bold", size = 10),
        panel.grid.minor = element_blank())

plot_dat <- clariplot3(mod.2)

p1 <- ggplot(plot_dat[["mfx.df"]], aes(x = monthsoffc, y = mean, group = factor(left), linetype = factor(left))) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = "grey70", alpha = 0.3) +
  geom_line() +
  facet_wrap(~ factor(populist)) +
  scale_linetype_manual(values = c(1, 2), labels = c("Right/Center", "Left")) +
  labs(x = NULL,
       y = "Spread volatility") +
  ggtitle(label = "EMBI Spread",
          subtitle = bquote("t-statistic on" ~ italic("Left Gov.")  %*% italic("Populist") %*% italic("Months in Office") ~ "is" ~ .(plot_dat[["t_int"]]))) +
  ylim(0, NA) +
  p_theme +
  NULL

plot_dat <- clariplot3(mod.4)

p2 <- ggplot(plot_dat[["mfx.df"]], aes(x = monthsoffc, y = mean, group = factor(left), linetype = factor(left))) +
  geom_ribbon(aes(ymin = lb, ymax = ub), fill = "grey70", alpha = 0.3) +
  geom_line() +
  facet_wrap(~ factor(populist)) +
  scale_linetype_manual(values = c(1, 2), labels = c("Right/Center", "Left")) +
  labs(x = NULL,
       y = "Spread volatility") +
  ggtitle(label = "CDS Spread",
          subtitle = bquote("t-statistic on" ~ italic("Left Gov.") %*% italic("Populist") %*% italic("Months in Office") ~ "is" ~ .(plot_dat[["t_int"]]))) +
  ylim(0, NA) +
  p_theme +
  NULL

# Combine plots

legend <- gtable_filter(ggplot_gtable(ggplot_build(p1 + theme(legend.position = "bottom"))), "guide-box")
lheight <- sum(legend$height)
p <- arrangeGrob(p1 + theme(legend.position = "none"), p2 + theme(legend.position = "none"),
                 ncol = 1, bottom = textGrob("Months in office", gp = gpar(fontfamily = "Myriad Pro")))
theight <- unit(12, "points")
p <- arrangeGrob(p, legend, heights=unit.c(unit(1, "npc") - lheight*2, lheight*1.5), ncol = 1)
grid.arrange(p)


# Subsample of democracies ------------------------------------------------

# Create criteria for inclusion:
# 1) Minimum country polity in the period >= 6
# 2) Average country polity in the period >= 6
# 3) Country with polity >= 6 in at least 75% of period covered

spread.dat <- spread.dat %>%
  group_by(iso3c) %>%
  mutate(min_polity = min(polity2, na.rm = TRUE),
         avg_polity = mean(polity2, na.rm = TRUE),
         polity_75pct = ifelse(quantile(polity2, probs = .75, na.rm = TRUE) >= 6, 1, 0),
         latam = ifelse(ht_region == 2, 1, 0)) %>%
  as.data.frame()

with(spread.dat, table(latam, polity_75pct))


fe.var <- "country" # Fixed-effect unit

y.var <- "d_embispread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.1 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$polity2 >= 6, ])
mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$polity2 < 6, ])

mod.3 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$polity_75pct == 1, ])
mod.4 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$polity_75pct == 0, ])

mod.5 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$avg_polity >= 6, ])
mod.6 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$avg_polity < 6, ])

mod.7 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$min_polity >= 6, ])
mod.8 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$min_polity < 6, ])


stargazer(mod.1[["cl.SE.table"]], mod.2[["cl.SE.table"]], mod.3[["cl.SE.table"]],
          mod.4[["cl.SE.table"]], mod.5[["cl.SE.table"]], mod.6[["cl.SE.table"]],
          mod.7[["cl.SE.table"]], mod.8[["cl.SE.table"]],
          add.lines = list(c("Observations",
                             mod.1[["nobs"]], mod.2[["nobs"]], mod.3[["nobs"]],
                             mod.4[["nobs"]], mod.5[["nobs"]], mod.6[["nobs"]],
                             mod.7[["nobs"]], mod.8[["nobs"]]),
                           c("Countries",
                             mod.1[["countries"]], mod.2[["countries"]], mod.3[["countries"]],
                             mod.4[["countries"]], mod.5[["countries"]], mod.6[["countries"]],
                             mod.7[["countries"]], mod.8[["countries"]]),
                           c("Country Fixed Effects",
                             "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("AIC",
                             mod.1[["aic"]], mod.2[["aic"]], mod.3[["aic"]],
                             mod.4[["aic"]], mod.5[["aic"]], mod.6[["aic"]],
                             mod.7[["aic"]], mod.8[["aic"]])),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)


# Regional heterogeneity --------------------------------------------------
# Is the effect specific to Latin America?

fe.var <- "country" # Fixed-effect unit

# Latin America

y.var <- "d_embispread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.1 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$ht_region == 2,])

y.var <- "d_cdsspread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$ht_region == 2,])

# Non-Latin America

y.var <- "d_embispread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.3 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$ht_region != 2,])

y.var <- "d_cdsspread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.4 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$ht_region != 2,])


stargazer(mod.1[["cl.SE.table"]], mod.2[["cl.SE.table"]], mod.3[["cl.SE.table"]], mod.4[["cl.SE.table"]],
          add.lines = list(c("Observations", mod.1[["nobs"]], mod.2[["nobs"]], mod.3[["nobs"]],
                             mod.4[["nobs"]]),
                           c("Countries", mod.1[["countries"]], mod.2[["countries"]], mod.3[["countries"]],
                             mod.4[["countries"]]),
                           c("Country Fixed Effects", "Yes", "Yes", "Yes", "Yes"),
                           c("AIC", mod.1[["aic"]], mod.2[["aic"]], mod.3[["aic"]],
                             mod.4[["aic"]])),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)

# Plot effect in Latin America vs Other Regions

# Function for plotting interaction in volatility
# equation using simulated Clarify CIs

clariplot4 <- function(mod){
  
  # Draw parameters from multivariate distribution
  S <- 10000
  set.seed(98765)
  sim.coef <- mvrnorm(n = S, mu = mod[["betas"]][,"coefficients"], Sigma = mod[["sandwich"]])
  sim.sigma.coef <- as.matrix(dplyr::select(as.data.frame(sim.coef), starts_with("sigma_")))
  
  # Set covariates at desired values
  monthsoffc <- seq(1, 96, by = 1)
  execleft.1 <- 1
  execleft.0 <- 0
  vix <- median(unique(spread.dat[, c("month_year", "d_vix")])$d_vix)
  kaopen <- median(spread.dat$kaopen, na.rm = TRUE)
  election_window <- 0
  
  out.1 <- matrix(NA, nrow = length(monthsoffc), ncol = 5)
  out.0 <- matrix(NA, nrow = length(monthsoffc), ncol = 5)
  
  for(i in 1:length(monthsoffc)){
    
    x.1 <- cbind(1, monthsoffc[i], execleft.1, vix, kaopen, election_window,
                 monthsoffc[i] * execleft.1)
    
    pred.sim.1 <- exp(sim.sigma.coef %*% t(x.1))
    
    x.0 <- cbind(1, monthsoffc[i], execleft.0, vix, kaopen, election_window,
                 monthsoffc[i] * execleft.0)
    
    pred.sim.0 <- exp(sim.sigma.coef %*% t(x.0))
    
    pred.sim <- pred.sim.1 - pred.sim.0
    
    out.1[i,] <- c(1, i, mean(pred.sim.1, na.rm = TRUE), quantile(pred.sim.1, c(0.025, 0.975), na.rm = TRUE))
    out.0[i,] <- c(0, i, mean(pred.sim.0, na.rm = TRUE), quantile(pred.sim.0, c(0.025, 0.975), na.rm = TRUE))
    
    rm(x.1, x.0, pred.sim.1, pred.sim.0, pred.sim)
  }
  
  out <- as.data.frame(rbind(out.1, out.0))
  colnames(out) <- c("left", "monthsoffc", "mean", "lb", "ub")
  
  # Get coefficient and t-stat on interaction term for annotation
  b_int <- signif(mod[["cl.SE.table"]]["sigma_monthsoffc:execleft", 1], digits = 2)
  t_int <- format(round(mod[["cl.SE.table"]]["sigma_monthsoffc:execleft", 3], digits = 2), nsmall = 2)
  
  return(list(out = out, b_int = b_int, t_int = t_int))
}

# Common plot theme

p_theme <- theme_bw(base_family = "Myriad Pro") +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        plot.title =  element_text(hjust = 0.5),
        plot.subtitle =  element_text(hjust = 0.5),
        #panel.grid.minor = element_blank(),
        panel.grid = element_line(linetype = 3))

plot_dat <- clariplot4(mod.1)

up_limit <- max(plot_dat[["out"]]$ub)

p1 <- ggplot(plot_dat[["out"]], aes(x = monthsoffc, y = mean, group = factor(left), linetype = factor(left))) +
  geom_ribbon(aes(x = monthsoffc, ymin = lb, ymax = ub), fill = "grey80", alpha = 0.3) +
  geom_line() +
  labs(x = NULL, y = NULL) +
  scale_linetype_manual(values = c(1, 2), labels = c("Right/Center", "Left")) +
  ylim(25, up_limit) +
  ggtitle(label = "Latin America",
          #subtitle = paste("Coef. on 'Left Gov. x Months in Office' is ", plot_dat[["b_int"]], "\n(t-statistic = ", plot_dat[["t_int"]], ")", sep = "")
          subtitle = bquote(atop("Coef. on" ~ italic("Left Gov.") %*% italic("Months in Office") ~ "is" ~ .(plot_dat[["b_int"]]), "t-statistic = " ~ .(plot_dat[["t_int"]])))) +
  p_theme+
  NULL

plot_dat <- clariplot4(mod.3)

p2 <- ggplot(plot_dat[["out"]], aes(x = monthsoffc, y = mean, group = factor(left), linetype = factor(left))) + 
  geom_ribbon(aes(x = monthsoffc, ymin = lb, ymax = ub), fill = "grey80", alpha = 0.3) +
  geom_line() +
  labs(x = NULL, y = NULL) +
  scale_linetype_manual(values = c(1, 2), labels = c("Right/Center", "Left")) +
  ylim(25, up_limit) +
  ggtitle(label = "Other Regions",
          #subtitle = paste("Coef. on 'Left Gov. x Months in Office' is ", plot_dat[["b_int"]], "\n(t-statistic = ", plot_dat[["t_int"]], ")", sep = ""))
          subtitle = bquote(atop("Coef. on" ~ italic("Left Gov.") %*% italic("Months in Office") ~ "is" ~ .(plot_dat[["b_int"]]), "t-statistic = " ~ .(plot_dat[["t_int"]])))) +
  p_theme +
  NULL

# Combine plots

legend <- gtable_filter(ggplot_gtable(ggplot_build(p1 + theme(legend.position = "bottom"))), "guide-box")
lheight <- sum(legend$height)
p <- arrangeGrob(p1 + theme(legend.position = "none"), p2 + theme(legend.position = "none"),
                 ncol = 2, left = "EMBI spread volatility", bottom = textGrob("Months in office", gp = gpar(fontfamily = "Myriad Pro")))
theight <- unit(12, "points")
p <- arrangeGrob(p, legend, heights=unit.c(unit(1, "npc") - lheight*2, lheight*1.5), ncol = 1)
grid.arrange(p)


# Baker & Greene 2011 party ideology scores for Latin America -------------

fe.var <- "country" # Fixed-effect unit

y.var <- "d_embispread" # Outcome measure

x.var <- c("monthsoffc")
x.int <- c("monthsoffc", "baker_ideol")
x.het <- c("monthsoffc")

mod.1 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "baker_ideol")
x.het <- c("d_vix", "kaopen", "election_window")

mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)


y.var <- "d_cdsspread" # Outcome measure

x.var <- c("monthsoffc")
x.int <- c("monthsoffc", "baker_ideol")
x.het <- c("monthsoffc")

mod.3 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "baker_ideol")
x.het <- c("d_vix", "kaopen", "election_window")

mod.4 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)


stargazer(mod.1[["cl.SE.table"]], mod.2[["cl.SE.table"]], mod.3[["cl.SE.table"]], mod.4[["cl.SE.table"]],
          add.lines = list(c("Observations", mod.1[["nobs"]], mod.2[["nobs"]], mod.3[["nobs"]], mod.4[["nobs"]]),
                           c("Countries", mod.1[["countries"]], mod.2[["countries"]], mod.3[["countries"]], mod.4[["countries"]]),
                           c("Country Fixed Effects", "Yes", "Yes", "Yes", "Yes"),
                           c("AIC", mod.1[["aic"]], mod.2[["aic"]], mod.3[["aic"]], mod.4[["aic"]])),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)


# Controlling for exchange rate movements ---------------------------------

fe.var <- "country" # Fixed-effect unit

y.var <- "d_embispread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation", "d_fxrate",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.1 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)


y.var <- "d_cdsspread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation", "d_fxrate",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)


stargazer(mod.1[["cl.SE.table"]], mod.2[["cl.SE.table"]],
          add.lines = list(c("Observations", mod.1[["nobs"]], mod.2[["nobs"]]),
                           c("Countries", mod.1[["countries"]], mod.2[["countries"]]),
                           c("Country Fixed Effects", "Yes", "Yes"),
                           c("AIC", mod.1[["aic"]], mod.2[["aic"]])),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)


# Controlling for sovereign credit ratings --------------------------------

fe.var <- "country" # Fixed-effect unit

y.var <- "d_embispread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation", "d_SPratingnumeric",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.1 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation", "d_moodysratingnumeric",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation", "d_fitchratingnumeric",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.3 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)


y.var <- "d_cdsspread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation", "d_SPratingnumeric",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.4 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation", "d_moodysratingnumeric",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.5 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation", "d_fitchratingnumeric",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.6 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)


stargazer(mod.1[["cl.SE.table"]], mod.2[["cl.SE.table"]], mod.3[["cl.SE.table"]],
          mod.4[["cl.SE.table"]], mod.5[["cl.SE.table"]], mod.6[["cl.SE.table"]],
          add.lines = list(c("Observations", mod.1[["nobs"]], mod.2[["nobs"]], mod.3[["nobs"]],
                                             mod.4[["nobs"]], mod.5[["nobs"]], mod.6[["nobs"]]),
                           c("Countries", mod.1[["countries"]], mod.2[["countries"]], mod.3[["countries"]],
                                          mod.4[["countries"]], mod.5[["countries"]], mod.6[["countries"]]),
                           c("Country Fixed Effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("AIC", mod.1[["aic"]], mod.2[["aic"]], mod.3[["aic"]],
                                    mod.4[["aic"]], mod.5[["aic"]], mod.6[["aic"]])),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)


# Controlling for political insitutions -----------------------------------

# Liberal democracy (V-Dem)
# Political constraints (Henisz)
# Central bank independence (Bodea & Hicks 2015)
# (Caution: slow-changing variables in fixed effects model)

fe.var <- "country" # Fixed-effect unit

y.var <- "d_embispread" # Outcome measure

x.var <- c("election_window", "v2x_libdem",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window", "v2x_libdem")

mod.1 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window", "polconiii",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window", "polconiii")

mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window", "lvaw",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window", "lvaw")

mod.3 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)


y.var <- "d_cdsspread" # Outcome measure

x.var <- c("election_window", "v2x_libdem",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window", "v2x_libdem")

mod.4 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window", "polconiii",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window", "polconiii")

mod.5 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window", "lvaw",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window", "lvaw")

mod.6 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)


stargazer(mod.1[["cl.SE.table"]], mod.2[["cl.SE.table"]], mod.3[["cl.SE.table"]], mod.4[["cl.SE.table"]],
          mod.5[["cl.SE.table"]], mod.6[["cl.SE.table"]],
          add.lines = list(c("Observations", mod.1[["nobs"]], mod.2[["nobs"]], mod.3[["nobs"]], mod.4[["nobs"]], mod.5[["nobs"]], mod.6[["nobs"]]),
                           c("Countries", mod.1[["countries"]], mod.2[["countries"]], mod.3[["countries"]], mod.4[["countries"]], mod.5[["countries"]], mod.6[["countries"]]),
                           c("Country Fixed Effects", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("AIC", mod.1[["aic"]], mod.2[["aic"]], mod.3[["aic"]], mod.4[["aic"]], mod.5[["aic"]], mod.6[["aic"]])),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)


# Comparing Chinn-Ito and Karcher-Steinberg KA openness -------------------

# Chinn-Ito

fe.var <- "country" # Fixed-effect unit

y.var <- "d_embispread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.1 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

y.var <- "d_cdsspread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

# Karcher-Steinberg

y.var <- "d_embispread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "ckaopen2010", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "ckaopen2010", "election_window")

mod.3 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

y.var <- "d_cdsspread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "ckaopen2010", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "ckaopen2010", "election_window")

mod.4 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

stargazer(mod.1[["cl.SE.table"]], mod.2[["cl.SE.table"]], mod.3[["cl.SE.table"]], mod.4[["cl.SE.table"]],
          add.lines = list(c("Observations", mod.1[["nobs"]], mod.2[["nobs"]], mod.3[["nobs"]], mod.4[["nobs"]]),
                           c("Countries", mod.1[["countries"]], mod.2[["countries"]], mod.3[["countries"]], mod.4[["countries"]]),
                           c("Country Fixed Effects", "Yes", "Yes", "Yes", "Yes"),
                           c("AIC", mod.1[["aic"]], mod.2[["aic"]], mod.3[["aic"]], mod.4[["aic"]])),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)


# Quarterly country-level macroeconomic indicators ------------------------

fe.var <- "country" # Fixed-effect unit

y.var <- "d_embispread" # Outcome measure

x.var <- c("election_window",
           "growth_qtr", "cab_qtr_gdp", "extdebt_gdp_qtr", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.1 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)


y.var <- "d_cdsspread" # Outcome measure

x.var <- c("election_window",
           "growth_qtr", "cab_qtr_gdp", "extdebt_gdp_qtr", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

stargazer(mod.1[["cl.SE.table"]], mod.2[["cl.SE.table"]],
          add.lines = list(c("Observations", mod.1[["nobs"]], mod.2[["nobs"]]),
                           c("Countries", mod.1[["countries"]], mod.2[["countries"]]),
                           c("Country Fixed Effects", "Yes", "Yes"),
                           c("AIC", mod.1[["aic"]], mod.2[["aic"]])),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)


# Month-year fixed effects ------------------------------------------------

fe.var <- "country" # Fixed-effect unit
time.var <- "month_year"

y.var <- "d_embispread" # Outcome measure

x.var <- c("monthsoffc")
x.int <- c("monthsoffc", "execleft")
x.het <- c("monthsoffc")

mod.1 <- hetreg(y.var, x.var, x.int, x.het, fe.var, time.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation", "kaopen") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, time.var, clusterSE = TRUE, data = spread.dat)


y.var <- "d_cdsspread" # Outcome measure

x.var <- c("monthsoffc")
x.int <- c("monthsoffc", "execleft")
x.het <- c("monthsoffc")

mod.3 <- hetreg(y.var, x.var, x.int, x.het, fe.var, time.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation", "kaopen") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.4 <- hetreg(y.var, x.var, x.int, x.het, fe.var, time.var, clusterSE = TRUE, data = spread.dat)

stargazer(mod.1[["cl.SE.table"]], mod.2[["cl.SE.table"]], mod.3[["cl.SE.table"]], mod.4[["cl.SE.table"]],
          add.lines = list(c("Observations", mod.1[["nobs"]], mod.2[["nobs"]], mod.3[["nobs"]], mod.4[["nobs"]]),
                           c("Countries", mod.1[["countries"]], mod.2[["countries"]], mod.3[["countries"]], mod.4[["countries"]]),
                           c("Country Fixed Effects", "Yes", "Yes", "Yes", "Yes"),
                           c("Month-Year Fixed Effects", "Yes", "Yes", "Yes", "Yes"),
                           c("AIC", mod.1[["aic"]], mod.2[["aic"]], mod.3[["aic"]], mod.4[["aic"]])),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)


# Latin America only w/ month-year FE -------------------------------------

fe.var <- "country" # Fixed-effect unit
time.var <- "month_year"

y.var <- "d_embispread" # Outcome measure

x.var <- c("monthsoffc")
x.int <- c("monthsoffc", "baker_ideol")
x.het <- c("monthsoffc")

mod.1 <- hetreg(y.var, x.var, x.int, x.het, fe.var, time.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation", "kaopen") 
x.int <- c("monthsoffc", "baker_ideol")
x.het <- c("d_vix", "kaopen", "election_window")

mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, time.var, clusterSE = TRUE, data = spread.dat)


y.var <- "d_cdsspread" # Outcome measure

x.var <- c("monthsoffc")
x.int <- c("monthsoffc", "baker_ideol")
x.het <- c("monthsoffc")

mod.3 <- hetreg(y.var, x.var, x.int, x.het, fe.var, time.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation", "kaopen") 
x.int <- c("monthsoffc", "baker_ideol")
x.het <- c("d_vix", "kaopen", "election_window")

mod.4 <- hetreg(y.var, x.var, x.int, x.het, fe.var, time.var, clusterSE = TRUE, data = spread.dat)

stargazer(mod.1[["cl.SE.table"]], mod.2[["cl.SE.table"]], mod.3[["cl.SE.table"]], mod.4[["cl.SE.table"]],
          add.lines = list(c("Observations", mod.1[["nobs"]], mod.2[["nobs"]], mod.3[["nobs"]], mod.4[["nobs"]]),
                           c("Countries", mod.1[["countries"]], mod.2[["countries"]], mod.3[["countries"]], mod.4[["countries"]]),
                           c("Country Fixed Effects", "Yes", "Yes", "Yes", "Yes"),
                           c("AIC", mod.1[["aic"]], mod.2[["aic"]], mod.3[["aic"]], mod.4[["aic"]])),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)


# Controlling for crises -----------------------------------------------------
# Laeven and Valencia crisis indicators
# Cruces and Trebesch debt restructuring indicator
# (EMBI data only; not enough data to estimate CDS models)

fe.var <- "country" # Fixed-effect unit

y.var <- "d_embispread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault",
           "sovereigndefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window", "sovereigndefault")

mod.1 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)


x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault",
           "sovereigndebtrestructure") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window", "sovereigndebtrestructure")

mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)


x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault",
           "currencycrisisstart") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window", "currencycrisisstart")

mod.3 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)


x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault",
           "sovereigndefault", "sovereigndebtrestructure", "currencycrisisstart") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window", "sovereigndefault", "sovereigndebtrestructure", "currencycrisisstart")

mod.4 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

stargazer(mod.1[["cl.SE.table"]], mod.2[["cl.SE.table"]], mod.3[["cl.SE.table"]], mod.4[["cl.SE.table"]],
          add.lines = list(c("Observations", mod.1[["nobs"]], mod.2[["nobs"]], mod.3[["nobs"]], mod.4[["nobs"]]),
                           c("Countries", mod.1[["countries"]], mod.2[["countries"]], mod.3[["countries"]], mod.4[["countries"]]),
                           c("Country Fixed Effects", "Yes", "Yes", "Yes", "Yes"),
                           c("AIC", mod.1[["aic"]], mod.2[["aic"]], mod.3[["aic"]], mod.4[["aic"]])),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)


# Temporal heterogeneity --------------------------------------------------

# Is the relationship different pre/post 2008 global financial crisis?

fe.var <- "country" # Fixed-effect unit

y.var <- "d_embispread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.1 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$year <= 2008,])

mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$year > 2008,])

stargazer(mod.1[["cl.SE.table"]], mod.2[["cl.SE.table"]],
          add.lines = list(c("Observations", mod.1[["nobs"]], mod.2[["nobs"]]),
                           c("Countries", mod.1[["countries"]], mod.2[["countries"]]),
                           c("Country Fixed Effects", "Yes", "Yes"),
                           c("AIC", mod.1[["aic"]], mod.2[["aic"]])),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)


# Controlling for close elections -----------------------------------------

# Create election window X close election interaction
spread.dat$election_windowXclose5 <- spread.dat$election_window * spread.dat$close_election5
spread.dat$election_windowXclose10 <- spread.dat$election_window * spread.dat$close_election10

fe.var <- "country" # Fixed-effect unit

y.var <- "d_embispread" # Outcome measure

x.var <- c("election_window", "close_election5", "election_windowXclose5",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window", "close_election5", "election_windowXclose5")

mod.1 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window", "close_election10", "election_windowXclose10",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window", "close_election10", "election_windowXclose10")

mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)


y.var <- "d_cdsspread" # Outcome measure

x.var <- c("election_window", "close_election5", "election_windowXclose5",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window", "close_election5", "election_windowXclose5")

mod.3 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

x.var <- c("election_window", "close_election10", "election_windowXclose10",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window", "close_election10", "election_windowXclose10")

mod.4 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat)

stargazer(mod.1[["cl.SE.table"]], mod.2[["cl.SE.table"]], mod.3[["cl.SE.table"]], mod.4[["cl.SE.table"]],
          add.lines = list(c("Observations", mod.1[["nobs"]], mod.2[["nobs"]], mod.3[["nobs"]], mod.4[["nobs"]]),
                           c("Countries", mod.1[["countries"]], mod.2[["countries"]], mod.3[["countries"]], mod.4[["countries"]]),
                           c("Country Fixed Effects", "Yes", "Yes", "Yes", "Yes"),
                           c("AIC", mod.1[["aic"]], mod.2[["aic"]], mod.3[["aic"]], mod.4[["aic"]])),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)


# Presidential vs parliamentary -------------------------------------------

# Is the effect specific to presidential systems?

spread.dat$presidential <- ifelse(spread.dat$system %in% c(0), 1, 0)
spread.dat$parliamentary <- ifelse(spread.dat$system == 2, 1, 0)

fe.var <- "country" # Fixed-effect unit

# Presidential

y.var <- "d_embispread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.1 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$presidential == 1,])

y.var <- "d_cdsspread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.2 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$presidential == 1,])

# Parliamentary

y.var <- "d_embispread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_embispreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.3 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$parliamentary == 1,])

y.var <- "d_cdsspread" # Outcome measure

x.var <- c("election_window",
           "cab_annual_gdp", "extdebt_gni_annual", "shortdebt_reserves_annual", "growth_annual", "d_inflation",
           "d_treasury10y", "d_vix", "d_energy_index", "d_eqprem", "d_cdsspreadREGIONdiff", "kaopen", "globaldefault") 
x.int <- c("monthsoffc", "execleft")
x.het <- c("d_vix", "kaopen", "election_window")

mod.4 <- hetreg(y.var, x.var, x.int, x.het, fe.var, clusterSE = TRUE, data = spread.dat[spread.dat$parliamentary == 1,])

stargazer(mod.1[["cl.SE.table"]], mod.2[["cl.SE.table"]], mod.3[["cl.SE.table"]], mod.4[["cl.SE.table"]],
          add.lines = list(c("Observations", mod.1[["nobs"]], mod.2[["nobs"]], mod.3[["nobs"]],
                             mod.4[["nobs"]]),
                           c("Countries", mod.1[["countries"]], mod.2[["countries"]], mod.3[["countries"]],
                             mod.4[["countries"]]),
                           c("Country Fixed Effects", "Yes", "Yes", "Yes", "Yes"),
                           c("AIC", mod.1[["aic"]], mod.2[["aic"]], mod.3[["aic"]],
                             mod.4[["aic"]])),
          dep.var.labels = "Sovereign Spread",
          no.space = TRUE)
